In [1]:
import plotly
import plotly.graph_objs as go
import plotly.express as px
from plotly.subplots import make_subplots

import torch
import numpy as np
import pandas as pd

import dill as pickle

import seaborn as sns
import random
import math
In [ ]:

Wave equation with one spatial variable

$$ \frac{\partial^{2} u}{\partial t^{2}} - \frac{1}{25} \frac{\partial^{2} u}{\partial x^{2}} = 0, $$ $$ \\ 71\times71, x \in [0; 1], t \in [0; 1]. $$

In [2]:
grid_res = 70
title = 'wave_equation'
df = pd.read_csv('C:\\Users\\Ksenia\\NSS\\ODE_projects\\SDE-DS\\wave_sln_70.csv', header=None)
initial_data = torch.from_numpy(df.to_numpy()).reshape(-1)

x_grid = np.linspace(0, 1, grid_res + 1)
t_grid = np.linspace(0, 1, grid_res + 1)

params = [x_grid, t_grid]

x = torch.from_numpy(x_grid)
t = torch.from_numpy(t_grid)

grid = torch.cartesian_prod(x, t).float()
In [3]:
u_sde = torch.load('C:\\Users\\Ksenia\\NSS\\ODE_projects\\SDE-DS\\average_wave_pt\\sols_stacked_35_solutions_shape_torch.Size([5041, 1]).pt', weights_only=True)
In [4]:
u_sde_mean_tens = u_sde.reshape(*[len(i) for i in params]).detach().cpu().numpy().T
In [5]:
u_main = torch.load('C:\\Users\\Ksenia\\NSS\\ODE_projects\\SDE-DS\\file_u_main_[35, 71, 71]_autograd_[0].pt', weights_only=False)
In [6]:
u_main.shape
Out[6]:
(35, 71, 71)
In [7]:
def calculate_statistics(u_main):

    mean_arr = np.zeros((u_main.shape[1], u_main.shape[2]))
    var_arr = np.zeros((u_main.shape[1], u_main.shape[2]))
    s_g_arr = np.zeros((u_main.shape[1], u_main.shape[2])) # population standard deviation of data.
    s_arr = np.zeros((u_main.shape[1], u_main.shape[2])) # sample standard deviation of data
    
    for i in range(u_main.shape[1]):
        for j in range(u_main.shape[2]):
            mean_arr[i, j] = np.mean(u_main[:, i, j])
            var_arr[i, j] = np.var(u_main[:, i, j])
            m = np.mean(u_main[:, i, j])
            s_arr[i, j] = math.sqrt(np.sum(list(map(lambda x: (x - m)**2, u_main[:, i, j])))/(len(u_main[:, i, j]) - 1))
    
    mean_tens = torch.from_numpy(mean_arr)
    var_tens = torch.from_numpy(var_arr)
    s_g_arr = torch.from_numpy(var_arr) ** (1/2)
    s_arr = torch.from_numpy(s_arr)
    
    # Confidence region for the mean
    upper_bound = mean_tens + 1.96 * s_arr / math.sqrt(len(u_main))
    lower_bound = mean_tens - 1.96 * s_arr / math.sqrt(len(u_main))
    
    mean_tens = mean_tens.reshape(-1)
    upper_bound = upper_bound.reshape(-1)
    lower_bound = lower_bound.reshape(-1)
    
    return mean_tens, upper_bound, lower_bound
In [8]:
u_bs_mean_tens, u_bs_upper_bound, u_bs_lower_bound = calculate_statistics(u_main)

Метод Соболя основан на разложении дисперсии функции $f(\mathbf{x})$ на вклады различных входных параметров. Формула выглядит следующим образом:

$$ f(\mathbf{x}) = f_0 + \sum_{i=1}^n f_i(x_i) + \sum_{1 \leq i < j \leq n} f_{ij}(x_i, x_j) + \cdots + f_{1,2,\ldots,n}(x_1, x_2, \ldots, x_n), $$ где:

SDE: Представляет компоненты, связанные с дисперсией первого порядка, т.е. с влиянием индивидуальных стохастических параметров.

BS: Преимущественно соответствует дисперсиям второго порядка, отражая условные вероятности и зависимости между параметрами.

$$ M[BS−M[SDE]] $$ BS — результат применения подхода на основе байесовской сети.

SDE — результат применения подхода на основе стохастических дифференциальных уравнений.

In [9]:
u_temp = u_main - u_sde_mean_tens

result_1 = np.mean(u_temp , axis=0)

print(result_1)

result_1.shape

[[-0.01292047 -0.01334606 -0.01373914 ...  0.01158872  0.0090557
   0.0050815 ]
 [ 0.04311887  0.0433047   0.04350033 ... -0.16430853 -0.1725123
  -0.18118833]
 [ 0.08509692  0.085596    0.08606815 ... -0.29792038 -0.31092852
  -0.32398838]
 ...
 [ 0.14299496  0.14392528  0.1447486  ... -0.5721952  -0.5990019
  -0.6262118 ]
 [ 0.07228713  0.07265554  0.07298737 ... -0.29356402 -0.31186715
  -0.3312171 ]
 [-0.01658184 -0.01707336 -0.01756203 ...  0.03265116  0.0267742
   0.01848607]]
Out[9]:
(71, 71)
In [10]:
# building 3-dimensional graph
fig = go.Figure(data=[
    go.Mesh3d(x=grid[:, 0], y=grid[:, 1], z=initial_data, name='Initial field',
              legendgroup='i', showlegend=True, color='rgb(139,224,164)',
              opacity=0.5),
    go.Mesh3d(x=grid[:, 0], y=grid[:, 1], z=result_1.reshape(-1), name='M[BS−M[SDE]]',
              legendgroup='s', showlegend=True, color='lightpink',
              opacity=1),
    go.Mesh3d(x=grid[:, 1], y=grid[:, 1], z=u_sde_mean_tens.reshape(-1), name='Solution field - M[SDE]',
              legendgroup='sde', showlegend=True, color='rgb(139,200,164)',),
])

fig.update_layout(scene_aspectmode='auto')
fig.update_layout(showlegend=True,
                  scene=dict(
                      xaxis_title='x1',
                      yaxis_title='x2',
                      zaxis_title='u',
                      zaxis=dict(nticks=10, dtick=1),
                      aspectratio={"x": 1, "y": 1, "z": 1}
                  ),
                  height=800, width=800
                  )

fig.show()

$$ M[BS]−M[SDE] $$

In [11]:
u_main_bs = np.mean(u_main, axis=0)
result_2 = u_main_bs - u_sde_mean_tens
print(result_2)
result_2.shape

[[-0.01292047 -0.01334606 -0.01373914 ...  0.01158872  0.0090557
   0.00508148]
 [ 0.04311887  0.04330469  0.04350033 ... -0.16430855 -0.17251235
  -0.18118832]
 [ 0.08509695  0.08559598  0.08606817 ... -0.29792035 -0.31092864
  -0.32398844]
 ...
 [ 0.14299497  0.1439253   0.14474854 ... -0.57219523 -0.5990019
  -0.62621176]
 [ 0.07228713  0.07265553  0.07298738 ... -0.29356405 -0.31186718
  -0.33121714]
 [-0.01658184 -0.01707336 -0.01756203 ...  0.03265116  0.02677419
   0.01848608]]
Out[11]:
(71, 71)
In [12]:
rmse = np.sqrt(np.mean((result_1 - result_2) ** 2))
print(f"Root Mean Squared Error (RMSE): {rmse}")

Root Mean Squared Error (RMSE): 3.5406208098720526e-07
In [13]:
# building 3-dimensional graph
fig = go.Figure(data=[
    go.Mesh3d(x=grid[:, 0], y=grid[:, 1], z=initial_data, name='Initial field',
              legendgroup='i', showlegend=True, color='rgb(139,224,164)',
              opacity=0.5),
    go.Mesh3d(x=grid[:, 0], y=grid[:, 1], z=result_1.reshape(-1), name='M[BS−M[SDE]]',
              legendgroup='s', showlegend=True, color='lightpink',
              opacity=1),
    go.Mesh3d(x=grid[:, 0], y=grid[:, 1], z=result_2.reshape(-1), name='M[BS]−M[SDE]',
              legendgroup='c', showlegend=True, color='blue',
              opacity=1),
    go.Mesh3d(x=grid[:, 1], y=grid[:, 1], z=u_sde_mean_tens.reshape(-1), name='Solution field - M[SDE]',
              legendgroup='sde', showlegend=True, color='rgb(139,200,164)',),
])

fig.update_layout(scene_aspectmode='auto')
fig.update_layout(showlegend=True,
                  scene=dict(
                      xaxis_title='x1',
                      yaxis_title='x2',
                      zaxis_title='u',
                      zaxis=dict(nticks=10, dtick=1),
                      aspectratio={"x": 1, "y": 1, "z": 1}
                  ),
                  height=800, width=800
                  )

fig.show()
In [20]:
max_value = np.max(result_1)
variance = np.var(u_sde_mean_tens)

print(f"max M[BS - M[SDE]]: {max_value}")
print(f"var SDE: {variance}")

max M[BS - M[SDE]]: 0.4558561146259308
var SDE: 3.770282030105591
In [18]:
x = np.linspace(-5, 5, 100)  
y = np.linspace(-5, 5, 100)  
fig = make_subplots(rows=1, cols=2, subplot_titles=('<b>Solution<b>', '<b>M[SDE]<b>'))
fig.add_trace(
    go.Heatmap(
        z=u_sde_mean_tens, x=x, y=y, colorscale='magma', opacity=1,
        showscale=True,
        colorbar=dict(x=1.0) 
    ),
    row=1, col=1
)
fig.add_trace(
    go.Heatmap(
        z=u_sde_mean_tens, x=x, y=y, colorscale='magma', opacity=1,
        showscale=False   
    ),
    row=1, col=2
)
font=dict(
        family="ATimes New Roman, Times, serif",  
        size=16,  
        color="black"  
    ),
title=dict(
        font=dict(
            family="Times New Roman, Times, serif",
            size=20,
            color="black"
         )
)
fig.show()
In [19]:
fig = make_subplots(
    rows=1, cols=1,
    subplot_titles=('<b>M[BN]−M[SDE]</b>'))

fig = go.Figure(
    data=go.Heatmap(
        z=result_2,
        x=x,
        y=y,
        colorscale='magma',
        opacity=1,
        showscale=True,
        colorbar=dict(x=1.02)
    )
)

fig.update_layout(
    title={
        'text': "<b>M[BN]−M[SDE]</b>",  
        'y': 0.9,                     
        'x': 0.5,                      
        'xanchor': 'center',           
        'yanchor': 'top',              },
    height=600,  
    width=600,  
    showlegend=False,
    xaxis=dict(scaleanchor="y"),  
    yaxis=dict(scaleanchor="x"),  
    xaxis2=dict(scaleanchor="y2"),  
    yaxis2=dict(scaleanchor="x2")
)
font=dict(
        family="ATimes New Roman, Times, serif",  
        size=16,  
        color="black" 
    ),
title=dict(
        font=dict(
            family="Times New Roman, Times, serif",
            size=20,
            color="black"
         )
)
fig.show()